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Abstract: We perform a global fit to the structure function F2 measured in lepton- 
proton experiments at small values of Bjorken-x, x < 0.01, for all experimentally 
available values of Q^, 0.045 GeV^ < < 800 GeV^ We show that the recent 
improvements resulting from the inclusion of running coupling corrections allow for a 
description of data in terms of non-linear QCD evolution equations. In this approach 
F2 is calculated within the dipole model with all Bjorken-x dependence described by 
the running coupling Balitsky-Kovchegov equation. Two different initial conditions 
for the evolution are used, both yielding good fits to data with x^/d.o.f. < 1.1. 
Data for the proton longitudinal structure function F^, not included in the fits, are 
also well described. Our analysis allows to perform a first principle extrapolation 
of the proton-dipole scattering amplitude. We provide predictions for F2 and Fl in 
the kinematical regions of interest for future colliders and ultra-high energy cosmic 
rays. A numerical implementation of our results down to x = 10~^^ is released as a 
computer code for public use. 
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1. Introduction 

The experimental data collected in electron-proton deep inelastic scattering (DIS) 
experiments [US] at small values of Bjorken-x constitute one of the most valu- 
able sources of information to test and explore the high-energy limit of Quantum 
Chromodynamics (QCD). The standard analyses (see [19] and references therein) 
of these data are usually made in the framework of fixed order DGLAP evolution 
equations in which various resummation schemes have also been essayed. On the 
other hand a description of available data in terms of the non-linear QCD evolution 
equations [20-27] has — despite phenomenological analyses (see e.g. [28-34]) being 
most suggestive of the presence of saturation effects, a crucial physical ingredient 
for the description of high-energy scattering in the small-x domain of DIS — been 
elusive so far. 

The saturation phenomenon is closely related to unitarity of the quantum field 
theory and is characteristic of dense partonic systems. It admits an intuitively clear 
physical picture in the infinite momentum frame. There, the gluon distribution func- 
tion xG{x, Q^) can be interpreted as the number of gluons in the proton wave function 
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localized within a transverse area inversely proportional to the photon virtuality Q^, 
and carrying a fraction of the proton longitudinal momentum x. For fixed Q^, the 
number of gluons in the proton wave function increases with decreasing x due to 
additional gluon emission or gluon branching. Such growth of gluon densities has 
been experimentally observed at HERA and, if extrapolated towards smaller values 
of X, would threaten the unitarity of the theory. Hence, the proton gets denser and 
gluon-gluon recombination processes, which are essentially non-linear, slow down 
the non-abelian avalanche towards small-x. This mechanism tames the subsequent 
growth of gluon densities, i.e. they saturate, thus preventing unitarity violations. 
The intrinsic momentum scale that determines the separation between the dilute 
and dense domains in the proton wave function is the saturation scale Q'^{x). This 
scale can be understood as the inverse transverse area inside which the probability of 
finding more than one gluon is of order one. It is a dynamic scale whose growth is de- 
termined by the interplay between the linear, radiative processes and the non-linear, 
recombination ones. 

All these qualitative ideas are cast in a definite theoretical framework, the Color 
Glass Condensate (CGC). The CGC is endowed with a set of pcrturbative, non- 
linear evolution equations, the Jalilian-Marian-Iancu-McLcrran-Wcigcrt-Leonidov- 
Kovner (JIMWLK) equation [20-25] and the Balitsky-Kovchegov (BK) equation 
[26,27], that describe the small- a; evolution of hadronic wave functions. However, 
rather than in terms of partonic densities, high-energy QCD evolution is more natu- 
rally formulated in terms of correlators of Wilson fines as effective degrees of freedom. 
The JIMWLK equation is equivalent to an infinite set of coupled, non-finear evo- 
lution equations for all correlators of the Wilson lines -also known as Balitsky's 
hierarchy. In the limit of large number of colors (Nc) the hierarchy reduces to a 
single equation - the BK equation- for the correlator of two Wilson lines or, equiv- 
alently, for the (imaginary part of the) dipole scattering amplitude A/". As we shall 
explain in detail in Section 2, in the dipole model the small- a; dependence of the 
different DIS cross sections is completely encoded in the dipole scattering amplitude, 
and thus describable by the JIMWLK-BK equations. 

Even though the JIMWLK equation comprises a richer physical input than the 
BK equation, the latter has become the most widely used tool to study the small- 
X dynamics. This is in part due to the relative simplicity of the BK equation with 
respect to JIMWLK, whose solution demands the use of rather complicated numerical 
methods [35]. Further, the difference between the solutions of the BK and JIMWLK 
equations turns out to be significantly smaller, of order 0.1% [35], than the a priori 
expected 0{1/N^) corrections. The origin of the smallness of the subleading-A'c 
corrections have been investigated recently in [36]. For these reasons, here we will 
consider the BK equation, rather than JIMWLK, as the starting point to analyse 
the experimental data on the proton structure functions at small-x. 

One of the first and most successful phenomenological applications of satura- 
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tion based ideas to the description of small-x DIS data is due to Golec-Biernat and 
Wiisthoff (GBW) [28]. Their pioneering work rehes on the use of the dipole model in 
QCD [37,38], together with a relatively simple model for the dipole-proton scattering 
amplitude encoding the basic features of saturation, to calculate the DIS total and 
diffractive lepton-proton cross sections. In particular, the proton saturation scale was 
parametrized as Q^{x) — {xq/x)^ GeV^. Fits to HERA data yielded Xq — 3 ■ 10~^ 
and A = 0.288. Several improvements of the GBW model for the dipole scattering 
amplitude were proposed later on in [29,31-34,39]. Very succinctly, some of these 
works [29] incorporated features of BFKL dynamics and explicit impact parameter 
dependence in the scattering amplitude [30,32,33], whereas [39] focused in including 
DGLAP evolution into the model, which resulted in a improved fit to the higher 
data. A first attempt of combining BK and DGLAP dynamics in the description of 
DIS data was made in [30] . Finally, the relation to heavy ion collisions was explored 
in [31,34]. Overall, these works reported an evolution speed compatible with the one 
obtained in the GBW model, A ~ 0.2 ^ 0.3. 

A natural question arises of why the BK-JIMWLK equations, the most solid 
theoretical tool available to describe the small-x dynamics of the dipole scattering 
amplitude and, in particular, the x-dependence of the saturation scale, have not been 
directly applied to the study of DIS small-x data. The answer to this question is 
given by the analytical [40,41] and numerical [42-44] studies of the leading-order 
(LO) BK equation. In these works the growth of the saturation scale yielded by 
the LO BK equation was determined to be ~ x~^^°, with Xlo — 4.88 A^cCKs/tt. 
Thus, the LO result predicts a much faster growth of the saturation scale (and hence 
of DIS structure functions) with decreasing x than the one extracted phenomeno- 
logically. This insufficiency of LO BK can only be circumvented by introducing an 
unreasonably small value for the fixed coupling, rendering any attempt to describe 
experimental data far from meaningful. 

It has been a long-standing expectation that higher order corrections to the orig- 
inal LO BK-JIMWLK equations could bring the theoretical predictions closer to ex- 
perimental observations. Indeed, numerical estimates for the running coupling [43,44] 
and energy conservation corrections [44,45] -both subleading physical contributions 
to the LO kernel- based on heuristic modifications of the LO kernel indicated a 
significant reduction of the evolution speed, thus pointing in the right direction. 
Moreover, running coupling effects appeared to dominate the contribution to the 
evolution kernel with respect to energy conservation effects [44]. However, it was 
not until recently that an explicit first principle calculation of the running coupling 
corrections to the evolution kernel was performed in [46-48] by including Nf cor- 
rections {Nf being the number of fiavors) into the evolution kernel to all orders and 
by then completing Nf to the one-loop QCD beta-function. The numerical study of 
the BK equation at all orders in agNf, performed in [49], reported a significant slow- 
down of the evolution speed with respect to the solutions of the LO equation, hence 
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rising the hopes that the improved equation might become a useful phenomenological 
tooL In its first successful application it was used to describe the energy and rapidity 
dependences of particle multiphcities produced in nucleus-nucleus collisions at the 
Relativistic Heavy Ion Collider (RHIC) at the BNL [50]. 

Significant progress has also been made recently in the determination of sublead- 
ing physical effects, other than running coupling corrections, to the LO BK equation, 
namely the inclusion of pomeron loops (see e.g. [51,52] and references therein), finite- 
Nc corrections [36] or the determination of the complete next-to-leading evolution 
kernel [53] to the BK equation. However, our current understanding indicates that 
the running coupling effects are dominant with respect to pomeron loops (or particle 
number fluctuations) [54] or finite corrections [36]. Wc will therefore limit our- 
selves in the present work to the analysis of DIS small-a; data via the BK equation 
including only running coupling corrections. 

The first goal of this paper is to prove the abihty of the BK equation including 
running coupling corrections to account for the small-x behavior of the total, F2, 
and longitudinal, Fl, structure functions measured in DIS experiments (a first step 
in this direction, yet unpublished, was reported in [55]). To that end we shall devise 
a global fit to the available experimental data with x < Xq = 10^^ and for all values 
of Q^. Analogously to previous works, our starting point will be the dipole model of 
QCD. The main novelty of our work is that the dipole-proton scattering amplitude, 
instead of being modeled, is calculated via numerical solutions of the BK equation 
including running coupling corrections. The free parameters in our fit, to be detailed 
in Section 2, are those related to the parametrization of the initial condition for the 
evolution, a global coefficient that sets the normalization and a constant which relates 
the running of the coupling in momentum space to that in dipole size. As we show in 
Section 4.1, the fits yield a good x^/d.o.f. < 1.1, thus demonstrating that such partial 
improvement of the LO BK equation suffices to reconcile the theoretical predictions 
with experimental results. Further, in Section 4.2 with all the free parameters fixed 
by the global fit of available data, we make predictions for the same observables 
at much smaller values of x. Such predictions are completely driven by non-linear 
QCD dynamics and could be directly tested at the proposed Electron-Ion CoUider 
(EIC) [56] or Large Hadron-electron CoUider (LHeC) [57] experimental facihties, 
where values of x as low as x ~ 10~^ for ~ 1 GeV^ could be reached. 

Second, the upcoming LHC experimental programs in proton-proton, proton- 
nucleus and nucleus-nucleus demand a detailed knowledge of hadronic wave functions 
or parton density functions (PDF) at very small input for the calculation 

of many different observables (see, for instance, the discussions in [19,58]). While 
global PDF fits provide a description of currently available data, additional the- 
oretical input is needed in order to safely extrapolate towards values of x so far 
unexplored empirically and for which additional saturation effects appear unavoid- 
able. A similar situation is found in cosmic rays physics [59, 60] , where the highest 
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center- of- mass energies reached in primary collisions are simply unattainable in ac- 
celerator experiments in the foreseeable future. In this work we set the ground for a 
systematic program oriented to provide parameter-free extrapolations of the dipole 
amplitudes (both for proton and nuclei) to very small values of x based on first 
principle calculations. Parametrizations of the dipole-proton scattering amplitudes 
down to very small x based on the results of this work are publicly available through 
simple numeric routines [61]. 



2. Setup 

In this section we briefly review, in a self contained manner, the main ingredients 
needed for the calculation of the inclusive and longitudinal DIS structure functions. 

2.1 Dipole model 

At X ^ 1, the inclusive structure function of DIS can be expressed as 

F2{x,Q^) = ^^—{aT + aL) , (2.1) 



em 



where ftem is the electromagnetic coupling and aT,L stands for the virtual photon- 
proton cross section for transverse (T) and longitudinal (L) polarizations of the 
virtual photon. The longitudinal structure function is obtained by considering only 
the longitudinal contribution: 

Fl{x, Q') = ar. . (2.2) 

It is well known that at high energies or small x (where the coherence length of the 
virtual photon fluctuation ^ {2mNx)^^ ~ 0.1/x fm ^ Rn, with itln and Rn 
the proton mass and radius respectively), and using light-cone perturbation theory, 
the total virtual photon-proton cross section can be written as the convolution of 
the light-cone wave function squared for a virtual photon to fluctuate into a quark- 
antiquark dipole, |^t,lP, and the imaginary part of the dipole-target scattering 
amplitude, Af. For transverse and longitudinal polarizations of the virtual photon 
one writes [37,38]: 

<JT,L(x,g') = ^ dz J dhdr\^T,L{^,Q^,r)\^Af{h,r,x), (2.3) 

where z is the fraction of longitudinal momentum of the photon carried by the 
quark, r is the transverse separation between the quark and the antiquark and b 
the impact parameter of the dipole-target collision (henceforth boldface notation 
indicates two-dimensional vectors). The wave functions |^r,Lp for the sphtting of 
the photon into a qq dipole are perturbatively computable within QED. We refer the 
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reader to e.g. [28] for explicit expressions to lowest order in ctem- All the information 
about the strong interactions - along with all x-dependence- in Eq. (2.3) is encoded 
in the dipole-proton scattering amplitude, J\f {h,r,x). Although this quantity is a 
genuinely non-perturbative object, its evolution towards smaller values of x can be 
studied perturbatively via the BK equation. On the contrary, its impact parameter 
dependence cannot be studied by means of the perturbative BK equation, since it is 
governed by long distance, non-perturbative physics. To circumvent this theoretical 
limitation we will resort to the translational invariance approximation (also used 
in [28]), which regards the proton as homogeneous in the transverse plane. Under 
this approximation the virtual photon-proton cross section Eq. (2.3) can be rewritten 
as follows: 

(Tt,l{^,Q^)^<Jo dz J dr|*T,L(-2,g',r)|W(r,y), (2.4) 

where r = |r| is the dipole size (the notation t> = |v| for all the 2-dimensional vectors 
will be also employed throughout the rest of the paper) and Uo is a dimensionful 
constant resulting from the b integration that sets the normalization - this will be 
one of the free parameters in our fits. Note that this result relies on the assumption 
that a factorized structure of x, r and b dependences remains unchanged through- 
out the evolution. In this case the parameter ao is related to the i-dependence in 
diffractive events, see e.g. [62]. On the other hand, this factorized structure may 
be assumed solely for the initial condition, while small-x evolution is performed, 
in the translational-invariant approximation, separately for every impact parameter 
(as done e.g. for nuclei in [63,64]). This results in a (Tq varying (increasing) with 
energy [65] . We leave this latter aspect for future studies. 

2.2 BK equation with running coupling 

The primary physical mechanism driving the small-a; evolution of the dipole scat- 
tering amplitude is the emission of soft gluons off either the quark or the antiquark 
in the original dipole. The leading order BK equation resumming the corresponding 
as ln(l/x) contributions to all orders reads 

dM{r, Y) 



dY 



J driK'^''{r,ri,r2) 
X [X{n,Y) +X(r2,Y)-X(r,Y) - Af{n,Y) X(r2,Y)] , (2.5) 

with the evolution kernel given by 

i^^°(r,ri,r,) = ^-f^, (2.6) 



27r2 



rr r. 



1 ' 2 



and r2 = r — Fi. Here, Y — \n{xo/x) is the rapidity variable and Xq is the value of x 
where the evolution starts, which should be small enough for the dipole model to be 
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applicable. In our case xq = 0.01 will be the highest experimental value of x included 
in the fit. 

The calculations in [46,47] proceeded by including OgNf corrections {Nf being 
the number of flavors) into the evolution kernel to all orders and by then completing 
Nf to the one-loop QCD beta-function via replacing Nf — 67r/52, with P2 = (H^c— 
2Nf)/{12'K). The calculation of the UgNf corrections is particularly simple in the 
s-channel light-cone perturbation theory (LCPT) formalism used to derive the BK 
and JIMWLK equations: there as Nf corrections are solely due to chains of quark 
bubbles placed onto the s-channel gluon lines, as sketched in Fig. lA. Importantly, 
at the same degree of accuracy a new physical channel is opened, namely the emission 
of a quark- ant iquark pair, instead of a gluon, as depicted in Fig. IB. The calculation 
in [48] relied instead on the use of dispersive methods, arriving at the same results 
as in the perturbative calculation in [47]. 

Neglecting the impact parameter dependence, the improved BK evolution equa- 
tion for the dipole scattering amplitude obtained after resumming the subleading 
UgNf corrections to all orders in [46, 47] can be written in the following, rather 
general form [49]: 

= n[M] - S[M] , (2.7) 

where both TZ and S are functionals of the dipole scattering amplitude, M . The 
first, running coupling, term 7^[A/'] in Eq. (2.7) gathers all the asNf factors needed 
to complete the QCD beta function, resulting in a functional form identical to the 
LO one but involving a modified kernel which provides the scale setting for the 
running of the coupling. In turn, the second term in the r.h.s. of Eq. (2.7), 5 [A/], 
the subtraction term, accounts for conformal, non running-coupling contributions. 



^ ^0 > ^0 




A 



Figure 1: Schematic representation of the diagrams contributing to the evolution to all 
orders in agNf. The s-channel gluon line can be attached to either the quark (upper line) 
or the antiquark (lower line). 

It would be erroneous to identify the gluon and quark-antiquark emission chan- 
nels with the running and subtraction terms in Eq. (2.7) respectively. Importantly, 
the quark-antiquark channel contains a logarithmic ultra-violet (UV) divergence as- 
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sociated to the emission of a zero size pair which, in the \aige-Nc hmit, is indistin- 
guishable from one gluon emission and therefore contributes to the running of the 
couphng on an equal footing. The emission of finite size quark-antiquark pairs is UV 
finite and does not contribute to the running of the coupling. Thus, the decompo- 
sition of the evolution kernel into running and subtraction contributions, although 
constrained by unitarity arguments, is not unique. This is due to the fact that there 
is some freedom in the way in which the UV divergence can be singled out from 
the conformal one, so in order to perform a decomposition like the one in Eq. (2.7) 
a precise separation scheme needs to be specified. Not surprisingly, the separation 
schemes proposed in [46] and [47] were different. For a detailed discussion on this 
subject we refer the reader to [49]. 

In this work we will consider only the running term in the evolution kernel. Ide- 
ally one would like to include the subtraction piece of the evolution kernel in practical 
applications as this would eliminate the uncertainty associated to the scheme choice 
and would provide a richer physical description of the small-x evolution of the dipole 
scattering amplitude. Unfortunately, its numerical evaluation [49] demands a very 
large computing time. For a global fit like the one presented in this work, in which 
the evolution is performed ~ 10^ times, such computing time is simply unaffordable. 
On the other hand, as shown in [49] the contribution to the complete evolution kernel 
stemming from the subtraction term is systematically smaller - and negligible at high 
rapidities - than the one arising from the running term. In particular, we will follow 
the prescription proposed by Balitsky in [46] to single out the running term since, as 
demonstrated in [49] , such choice minimizes the contribution to the evolution of the 
subtraction term, neglected in what follows, with respect to the separation scheme 
proposed in [47]. 

Finally, after dropping the subtraction term from Eq. (2.7), the BK evolution 
equation including only running coupling corrections reads 

^%^ = KB-M, (2.8) 

where the running coupling functional is identical to the LO equation: 

n^^\M] = j driir^'^'(r,ri,r2) 

X [jV(ri, Y) + M{r2, Y) - M{r, Y) - M{r^, Y) M{r2, Y)] , (2.9) 

but with a modified evolution kernel that includes running coupling corrections. 
Using Balitsky's prescription, the kernel for the running term reads [46] 

r2 1 faJr^) \ 1 /aJr^) \1 ^ 

^ + - - 1 + - -n ■ 2-10 

rfr| rf \as{ri) J r| \as{rf) J 



i^^'^Xr,ri,r2) 



NcaA 
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2.3 Regulcirization of the infrared dynamics 

The BK equation is an integro-differential equation that involves integration over 
all available phase-space for soft gluon emission. In the running coupling case, Eqs. 
(2.8-2.10), the coupling has to be evaluated at arbitrarily large values of the dipole 
size (small gluon momentum), and a regularization prescription to avoid the Landau 
pole becomes necessary. A celebrated feature of the BK equation is its ability to 
fix [66] the problem of infra-red diffusion characteristic of its linear counterpart, the 
BFKL equation. The non-linear terms in the BK equation ensure that the dynamics 
in the phase space region within the unitarity limit, i.e. for r ^ 1/Qs, is frozen. 
Such feature is shared by both the LO and running coupling BK equations, since it 
is ultimately rooted in the non-linear combination of AA's in the r.h.s. of Eq. (2.9), 
which is identical in both cases. Thus, if Qg is perturbatively large, Qs ^ ^qcd, 
then all the relevant dynamics takes place deep in the ultra-violet region of the phase 
space, r < 1/Qs- In such scenario the details about the regularization of the running 
coupling in the infra-red become irrelevant for the result of the evolution. 

Unfortunately, we can anticipate that such will not be the case in this work. 
Taking the results by Golec-Biernat and Wiisthoff [28] as a guidance, one can estimate 
that the proton saturation scale at the largest values of Bjorken-x to be considered in 
this work, X ~ 10-2, is of the order oi QHx^lO-"^) ^ {S-lQ-yiO-y'^^^ GeV^ ~ 0.36 
GeV^. The fits to be presented in Section 4 yield even smaller values of the initial 
saturation scale of the proton. Although larger than Ag^^, such values for the 
initial scale are not large enough to avoid sensitivity to the infra-red (IR) dynamics. 
Actually, the detailed study of the infrared-renormalon ambiguities carried out in [48] 
demonstrated that the sensitivity of the solutions of the evolution equation to several 
different prescriptions used to regularize the coupling is relatively large even for initial 
saturation scales as large as ~ l-i-2 GeV^. On the bright side, theoretical studies 
of the Schwinger-Dyson equations for the gluon propagator in the IR and lattice 
QCD results (see e.g. [67,68] and references therein) indicate that the strong coupling 
freezes to a constant value, a/r, in the IR. Moreover, the value at which the coupling 
freezes has been determined to he ajr ~ 0.5 -i- 0.7. While these results are somewhat 
controversial and yet subject to discussion in the literature, in particular the very 
definition of an infrared coupling, we will take them as a guidance to regularize the IR 
dynamics. Otherwise, our prescription can be regarded as purely phenomenological. 

Thus, for small dipole sizes r < r/^, with ^^(ry^) = afr — 0.7, we shall evaluate 
the running coupling according to the usual one-loop QCD expression: 



with Nf — 3, whereas for larger sizes, r > r/^, we freeze the coupling to the fixed value 
afr = 0.7. We take Aqcd = 0.241 GeV, such that as{Mz) = 0.1176, with Mz the 
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mass of the Z boson. The factor under the logarithm in Eq. (2.11) will be one of 
the free parameters in the fit. It refiects the uncertainty in the Fourier transform from 
momentum space, where the original calculation of agN f corrections was performed, 
to coordinate space. Alternatively, we could have fixed to the value suggested 
in [47], er^l'^-'^^iE ^ q^^A chosen either Kqcd or a^r as the free parameters controlling 
the IR dynamics. Indeed, we have checked that such choices yield equally good fits 
as those presented in Section 4 without changing much the value of the other free 
parameters. However, both a^r and, specially, ^qcd-, are more tightly constrained 
from both theoretical and phenomenological studies than C^. 

2.4 Initial conditions for the evolution 

Finally we have to specify the initial condition (i.e.) for the evolution or, equivalently, 
the precise shape of the proton unintegrated gluon distribution (UGD), (f){x,k), at 
the highest experimental value of Bjorken-x included in the fit, xq = 0.01 (which, 
by definition, corresponds to rapidity Y — 0). The UGD is related to the dipole 
scattering amplitude via a Fourier transform: 



This is a genuinely non-perturbative object which needs to be modeled. We will 
consider two different families of initial conditions. The first one is inspired in the 
original GBW ansatz [28] for the dipole scattering amplitude and parametrized in 
the following way: 



The second family of initial conditions [69] follows closely the McLerran-Venugopalan 
(MV) model: 



where Q% is the initial saturation scale squared in both cases. 

Eqs. (2.13) and (2.14) differ with respect to the original GBW and MV models 
in the inclusion of an anomalous dimension, 7, which will be another of the free 
parameters in the fit. The GBW and MV functional forms are recovered by setting 
7 = 1 in Eq. (2.13) and Eq. (2.14) respectively. The anomalous dimension controls 
the slope of the scattering amplitude in the transition from the dilute region to the 
black disk region. The main difference between MV and GBW i.e. is their different 
UV behaviour, which is more easily appreciated in momentum space. For 7 = 1 
and large transverse momenta k, the UGD resulting from the MV i.e. falls off as 
(P^y as expected from rather general perturbative considerations, while the 
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(2.13) 




(2.14) 
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GBW i.e. falls off exponentially, (jp^^ ~ exp {—k'^/Ql). It is well known that the 
solutions of the BK equation, both at LO and including higher order corrections, do 
not respect the relatively simple functional forms in Eq. (2.13) and Eq. (2.14). On 
the contrary, they can be roughly characterized by an r- and F-dependent anomalous 
dimension, 7(r, F), with 7 ^ 1 for r — > 0. Clearly a constant value of 7 7^ 1 would 
not respect such condition. However, the main contribution to the DIS cross section 
given by Eq. (2.4) originates from the region 1/Q < r < 1/Qs- The contribution 
from the dilute U V region r < l/Q is much smaller and therefore we will not consider 
additional refinements of the initial conditions in Eq. (2.13) and Eq. (2.14), which 
would come at the prize of adding new, spurious parameters into the fit. [Actually, 
the results of the fit shows that for the GBW i.e. the preferred value is 7 = 1, so it 
will be fixed for this initial condition.] Finally, the constant term under the logarithm 
in the MV initial condition, e, has been added to regularize the exponent for large 
values of r. 

2.5 Summciry of the theoretical setup and free peirameters 

In summary, we will calculate the total DIS inclusive and longitudinal structure func- 
tions according to the dipole model under the translational invariant approximation 
Eq. (2.4). The small- a; dependence is completely described by means of the BK equa- 
tion including running coupling corrections, Eqs. (2.8-2.10), for which two different 
initial conditions GBW and MV, Eqs. (2.13) and (2.14), are considered. All in all, 
the free parameters to be fitted to experimental data are: 

• (7o : The total normalization of the cross section in Eq. (2.4). 

• (5^0 ■ The saturation scale of the proton at the highest experimental value of 
Bjorken-x included in the fit, Xq = 10~^, in Eqs. (2.13) and (2.14). 

• C^: The parameter relating the running of the coupling in momentum space 
to the one in dipole size in Eq. (2.11). 

• 7 : The anomalous dimension of the initial condition for the evolution in Eqs. 
(2.13) and (2.14). As discussed in Section 4, in some cases (GBW) 7 can be 
fixed to 1, obtaining equally good fits to data than when it is considered a free 
parameter. 

3. Numerical method and experimental data 

The experimental data included in the fit to F2{x,Q^) have been collected by the 
E665 [1] (FNAL), the NMC [2] (CERN-SPS), the HI [3-8] (HERA) and the ZEUS 
[9-16] (HERA) experimental Collaborations. We have considered data for x < 10~^ 
and for all available values of Q'^, 0.045 GeV^ < < 800 GeV^. 
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The only published direct measurements of the longitudinal structure function 
Fl{x,Q^) were obtained recently by the HI [17] and ZEUS [18] Collaborations, and 
they are not included in the fit. 

All in all, 847 data points are included. Statistical and systematic uncertain- 
ties were added in quadrature, and normalization uncertainties not considered. [A 
more involved treatment separating uncorrelated and correlated/normalization er- 
rors could be done only at the expense of adding one more fitting parameter for 
each of the 17 data sets used, thus making the minimization task impossible due to 
CPU-time requirements.] Since the minimization algorithms require a large number 
of calls to the function we have implemented a parallelization of the numeric code. 
Finally, the BK evolution equation including running coupling corrections is solved 
using a Runge-Kutta method of second order with rapidity step Ahy = 0.05, see 
further details in [49]. 

In order to smoothly go to photoproduction, we follow [28] and use the redefini- 
tion of the Bjorken variable 

with ruf — O.IA GeV for the three hght flavors we consider in Eq. (2.4). 
4. Results 

4.1 Fits to F2 and description of 

The values of the free parameters obtained from the fits to data for the two different 
initial conditions, GBW and MV, are presented in Table 1. A partial comparison 
between the experimental data [1-16] and the results of the fit for F2{x, Q^) is shown 
in Fig. 2. 



Initial condition 


ctq (mb) 


Q'so (GeV^) 


C2 


7 


xVd.o.f. 


GBW 


31.59 


0.24 


5.3 


1 (fixed) 


916.3/844=1.086 


MV 


32.77 


0.15 


6.5 


1.13 


906.0/843=1.075 



Table 1: Values of the fitting parameters from the fit to F2{x,Q^) data from [1-16] with 
X < 10-2 and for all available values of Q^ 0.045 GeV^ < < 800 GeV^. 



On the other hand, Fl{x, Q"^) offers an additional constrain on the gluon distri- 
bution and is expected to have more discriminating power on different approaches, 
particularly in the \ow-Q^ region [70]. In Fig. 3 we show a comparison between 
experimental data [17, 18] and our predictions for Fl{x, Q^). 

Several comments are in order. First, the two different initial conditions yield 
very good fits to F2-data, with x^/d.o.f. ~ 1, and almost identical results for Fj,. 
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Figure 2: Comparison between a selection of experimental data [1-16] and the results 
from the fit for F2{x, Q^). Solid red lines correspond to GBW i.e., and dotted blue ones to 
MV i.e. The error bars correspond to statistical and systematic errors added in quadrature. 



As remarked in the previous Section the main difference between the two initial 
conditions is their behavior at small r. In principle this difference is large, but the 
fact that the values of 7 resulting from the fit are different for the different initial 
conditions, should compensate it in a limited region of r. We thus conclude that the 
kinematical coverage of the existing experimental data on F2 (and Fl) is too small 
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Figure 3: Comparison between experimental data from the HI [17] (upper plot) and 
ZEUS [18] (lower plot) Collaborations and the predictions of our model for Fl{x, Q^). Red 
solid lines and open squares correspond to GBW i.e., and blue dotted lines and open circles 
to MV i.e. The theoretical results have been computed at the same (x) as the experimental 
data, and then joined by straight lines. The error bars correspond to statistical and sys- 
tematic errors added in quadrature for those data coming from [17], while they correspond 
to the error quoted for the unconstrained fit for those data coming from [18]. 



to allow a discrimination of the different UV behaviors of the two employed i.e. 
Second, the fits using GBW i.e. and obtained by letting 7 vary as a free pa- 
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rameter, do not show an improvement with respect to those obtained by fixing it to 
7 = 1. On the contrary, the fits using MV i.e. do improve by letting 7 be a free 
parameter, which takes a value slightly larger than one, 7 = 1.13. 

Third, although the two different fits yield pretty different values of the initial 
proton saturation scale, this apparent discrepancy is due to the different functional 
forms for GBW and MV i.e. If we redefine the initial saturation scale for the MV i.e. 
via the condition N^^{r = l/(5',o,My> ^ = 0) = 1 - e"^/'' (see Section 4.3), we will 
get Q'sq^mv ~ O-l^ GeV^, which is closer to the GBW result. Therefore we conclude 
from our study that the saturation scale of the proton, obtained in our fit within the 
dipole model (considering only three active fiavors and translational invariant initial 
conditions i.e. a proton with a constant profile) at a; = 0.01, lies in the range 

0.19 GeV^ < QIq < 0.25 GeV^ 

Fourth, the values of cjq obtained from the fits are very close to each other. 
This supports the assumption of translational invariance. Furthermore, the obtained 
values of ctq — 32 mb correspond, assuming a Gaussian form factor for the proton 
[62], to a diffractive exponential slope cro/(47r) ~ 6.5 GeV~^ in agreement with 
experimental data [71], see the comments below Eq. (2.4). 

Fifth, we have checked that the quahty of the fit and the values of the parameters 
are stable under the restriction of the data range to the region < 50 GeV^ (which 
leaves 703 data points for the fit) . While in principle the dipole model should be more 
suitable for the description of structure functions in the region of low and moderate 
Q^, we take this stability as a signal that there is no tension in the fit with the 
large-Q^ data. 

Finally, the agreement of our predictions for Fl{x,Q'^) with the experimental 
data [17, 18] is of the same quality as other based on fixed-order NLO and NNLO 
DGLAP evolution, see the comparison in [17,18]. As discussed in [70], data at smaller 
may offer the possibility of discriminating different approaches. 

4.2 Predictions for future experimental programs 

Besides available experimental data, the experimental programmes at the LHC will 
test [19, 58, 72] our understanding of the small-x behavior of the nucleon structure. 
There are also proposals of future lepton-hadron colliders [56,57] in which new mea- 
surements of structure functions at smaller x would be performed. Further, the 
physics of high-energy cosmic rays is expected to be infiuenced by small-x phenom- 
ena [59,60]. Therefore, we find it worth to show in Fig. 4 our predictions for F2 and 
Fl in a broad, yet experimentally unexplored region of x and Q^. 

Two facts should be highlighted. First, the striking agreement of the predictions 
-which makes them more reliable- from both initial conditions. Second, that at 
large Q'^ and small x the effect of saturation, namely the flattening of the structure 
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Figure 4: Predictions for F2{x,Q'^) (top) and Fl{x,Q'^) (bottom) versus x, for 10 ^ < 
X < 10^2 and = 10^\ 1, 10, 10^, 10^ GeV^ (lines from bottom to top). Solid black lines 
show the results obtained with GBW i.e., and dotted red lines those obtained with MV i.e. 

function, is more apparent in Fl than in F2. This fact stresses, in our view, the 
importance of Fl measurements to distinguish different scenarios for the small-x 
dynamics: fixed order perturbative QCD, resummation schemes or saturation models 
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[70]. 



4.3 Parametrizations of the dipole-proton scattering amplitude 

With all the uncertainties associated to the initial condition for the evolution fixed by 
the fit to F2 presented in the previous sections, we can now evolve the proton-dipole 
scattering amplitude to much smaller values of x. Such extrapolation is completely 
driven by small-x evolution including running coupling corrections and can be used 
to calculate several different observables relevant for the LHC and cosmic ray physics. 
We have performed the evolution down to a; = 10~^^. The resulting proton-dipole 
scattering amplitude is plotted in Fig. 5 for three values of x {x = 10^^,5 ■ 10~^ 
and 5 ■ 10^^) both for MV and GBW i.e. and has been made public through simple 
fortran routines in [61]. From the solutions of the evolution in Fig. 5 we can extract 
the proton saturation scale Qs{x) through the condition 



It is important to note that the values of Qs{x) presented in Fig. 6 are dependent 
on the choice of k, in Eq. (4.1). Following the original GBW prescription we take 



Different choices of k can affect the numerical value of Qs{x) by a factor ~ 2 3. 
Keeping in mind such ambiguity in its extraction from the numerical solutions of the 
evolution equation, we can estimate the value of the proton saturation scale at LHC 
energies. Using 2 — > 1 kinematics to compute the smallest value of Bjorken-x probed 
in proton-proton collisions, x = (2M/y^)e~^, where M is the invariant mass of the 
produced system (one hadron, dileptons,...), ^/s = 14 TeV is the collision energy 
and y the rapidity of the produced particle, we get (fixing M = 1 GeV) that the 
saturation scale of the backward-moving proton at the LHC at rapidities y = 0,3 
and 6 is — 0.55-^-0.7, 1.3-^-1.7 and 3 -=-4 GeV^ respectively. Such values are large 
enough to suggest that saturation effects in proton-proton collisions at the LHC may 
be detectable, specially at forward rapidities. 

5. Conclusions 

We presented a new approach towards a systematic quantification of parton distri- 
butions at small-a; directly in terms of non-linear QCD evolution equations. This 
approach has become feasible thanks to the recent calculation of the running cou- 
pling corrections to the BK equation. In this work we performed a global fit to the 
available experimental data for F2{x, Q"^) measured in electron- proton scattering for 
x < 10~^ and all values of Q^. The calculation of the structure functions F2 and 
Fl is done within the dipole model under the translational invariant approximation 



J\f{r^l/Qs{x),x)^K^O{l). 



(4.1) 



K = l-exp[-l/4] ~0.22. 



(4.2) 
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Figure 5: Dipole scattering amplitude obtained from the fits for the two different initial 
conditions, MV (dashed blue) and GBW (solid red) at x = 10~^, 5- 10~^ and 5- 10~^ (from 
right to left). 

and considering just three active flavors. The main novehy of this work with respect 
to previous phenomenological analyses is the direct use of the running coupling BK 
equation to describe the small-x dependence of the structure functions. We find a 
very good agreement with experimental data with only three (four) free parameters 
using GBW (MV) initial conditions for the evolution. Available data on F^, not 
included in the fit, are also well described. We present predictions for both F2 and 
Fl in the kinematic regime relevant for future accelerators and ultra high-energy 
cosmic rays. We also provide the evolved proton-dipole scattering amplitude down 
to values of x = 10"^^ through a simple computer code for public use [61]. Further 
extension of this work to nuclear targets and hadronic and nuclear collisions is under 
way. 

In conclusion, we find that the recent progress in our knowledge of non-linear 
small-x evolution brings us to an unprecedented level of precision allowing for a direct 
comparison with experimental data. This provides a solid theoretical extrapolation 
of parton densities towards yet empirically unexplored kinematic regions. 
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